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Abstract: 

In many applications (hypergeometric-type) special functions like orthogonal polynomials are 
needed. For example in more than 50 % of the published solutions for the (application-oriented) 
questions in the "Problems Section" of SIAM Review special functions occur. 

In this article the Mathematica package SpecialFunctions which can be obtained from the URL 
http://www.zib.de/koepf is introduced [15]. Algorithms to convert between power series rep- 
resentations and their generating functions is the main topic of this package ([8]-[15]), extending 
the previous package PowerSeries [12]. Moreover the package automatically finds differential 
and recurrence equations ([13]-[14]) for expressions and for sums (the latter using Zeilberger's 
algorithm ([23], [18], [13]). 

As an application the fast computation of polynomial approximations of solutions of linear dif- 
ferential equations with polynomial coefficients is presented. This is the asymptotically fastest 
known algorithm for series computations, and it is much faster than Mathematica's builtin 
Series command if applicable. Many more applications are considered. 

Finally the package includes implementations supporting the efficient computation of classical 
continuous and discrete orthogonal polynomials. 

1 Holonomic Functions 

A homogeneous linear differential equation 

m 

with polynomial coefficients Pk{x) is called holonomic, as is the corresponding f(x). 
Holonomic functions have nice algebraic properties: In 1980 Stanley [22] proved by algebraic 
arguments that the sum and product of holonomic functions, and their composition with 
algebraic functions, in particular with rational functions and rational powers, form holonomic 
functions, again. 

By iterative differentiation and the use of Gaussian elimination, one can construct the resulting 
differential equations algorithmically. 

The funny thing is that these algorithms had been known already in the last century ([l]-[2]), 
but because of their complexity had fallen into oblivion. 

The package SpecialFunctions contains implementations of these algorithms. The Mathe- 
matica procedure SumDE [DEI , DE2 , f [x] ] computes the holonomic differential equation for the 
sum of two functions which correspond to the holonomic differential equations DEI and DE2, 
written in terms of f [x] . Here are some examples: 
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In[l]:= <<SpecialFunctions ' 

SpecialFunctions , (C) Wolfram Koepf , version 1.00, November 20, 1996 
Fast Zeilberger, (C) Peter Paule and Markus Schorn (V 2.2) loaded 



In [2] : = SumDE [f ' ' [x] +f [x] ==0 , f ' [x] -f [x] ==0 , f [x] ] 



(3) 



Out [2] = f [x] - f ' [x] + f " [x] - f 



[x] == 



In [3] : = ProductDE [f ' ' [x] +f [x] ==0 , f ' [x] -f [x] ==0 , f [x] ] 



Out [3] = 2 f [x] - 2 f ' [x] + f " [x] == 



In [4] : = SumDE [f " [x] -x*f [x] ==0 , f ' [x] -f [x] ==0 , f [x] ] 



2 

Out [4]= (1 - x + x ) f[x] + (1 - x) x f'[x] - 



(3) 



> 



x f " [x] + (-1 + x) f 



[x] == 



In [5] : = ProductDE [f ' ' [x] -x*f [x] ==0 , f ' [x] -f [x] ==0 , f [x] ] 
Out [5]= (-1 + x) f[x] + 2 f'[x] - f"[x] ==0 

In every day language the first computation reads as: The sum of the functions sin a; (satisfying 
f"{x) + f(x) = 0) and e x (satisfying f'(x) — f(x) = 0), sin a; + e x , satisfies the differential 
equation f(x) — f'{x) + f"(x) — f"'(x) = 0. Actually any linear combination a sin x + b cos x + 
c e x satisfies this resulting differential equation. In the next line the differential equation for 
the product sin^e^ (or cosrre^) is computed. Note that the further computations are rather 
similar, although one of the corresponding functions is much more advanced, namely an Airy 
function (satisfying f"(x) — x f(x) = 0). 

If functions are given by expressions, holonomic differential equations can be determined as 
well. The procedure HolonomicDE[expr,f [x]] internally uses the sum and product algo- 
rithms: 

In [6] : = HolonomicDE [AiryAi [x] *Exp [x] , f [x] ] 

Out [6]= (-1 + x) f[x] + 2 f'[x] - f"[x] ==0 

In [7]:= HolonomicDE [Exp [alpha x]*Sin[beta x],f[x]] 



Out [7] 



2 2 

(alpha + beta ) f [x] - 2 alpha f ' [x] + 



> 



f " [x] == 



In [8] := de=HolonomicDE[ArcSin[x] ,f [x]] 

2 

Out [8]= x f ' [x] + (-1 + x ) f " [x] == 



In [9] : = Nest [ProductDE [de , # , f [x] ] & , de , 4] 



2 



2 

Out [9]= x f'[x] + (-16 + 31 x ) f"[x] + 
2 (3) 

> 15 x (-5 + 6 x ) f [x] + 

2 (4) 

> 5 (1 - x) (1 + x) (4 - 13 x ) f [x] + 

2 2 (5) 

> 15 (-1 + x) x (1 + x) f [x] + 

3 3 (6) 

> (-1 + x) (1 + x) f [x] == 

The last statement computes the holonomic differential equation for arcsin 5 x. This result can 
also be obtained by the simple statement HolonomicDE[ArcSin[x] ~5,f [x]] . 

2 Series Representations 

By an easy procedure (equating coefficients) any holonomic differential equation can be con- 
verted to a holonomic recurrence equation for a series representation 

oo 

f(x) = a k xk ■ 

k=—oo 

It turns out that a function f(x) is holonomic if and only if its series coefficients are 
holonomic, i.e., if the latter satisfy a homogeneous, linear recurrence equation 

M 

^2 ak{n) f n +k{x) = 

with polynomial coefficients in n. The procedure DEtoRE[DE,f [x] ,a[k]] converts the 
differential equation DE for f(x) to the corresponding holonomic recurrence equation for the 
coefficients a&, and REtoDE[RE,a[k] , f [x]] does the corresponding back conversion. We have 
for example 

In [10] := de=HolonomicDE[ArcSin[x] ~2,f [x] ] 
Out [10]= f ' [x] + 3 x f " [x] + 

(3) 

> (-1 + x) (1 + x) f [x] == 
In[ll]:= re=DEtoRE[de,f [x] ,a[k]] 

3 

Out [11]= k a[k] - k (1 + k) (2 + k) a[2 + k] == 
In [12]:= REtoDE[re,a[k] ,f [x]] 



3 



(3) 2 (3) 

Out [12]= f ' [x] + 3 x f " [x] - f [x] + x f [x] == 

If for a given input function f(x) the resulting recurrence equation has only two summands, one 
has a generalized hypergeometric function (see e.g. [21]), and the coefficients can be computed 
explicitly. 

The generalized hypergeometric series is given by 



/ ai a 2 ••• % 
P q h b 2 ■■■ b q 



A;=0 



x k 



(ai)k ■ ( a 2)k ■ ■ ■ (a p )k % k 



k=0 



h)k ■ (h)k- ■ ■ (b q )kk\ 



k 

where (a)k ■= Yl {o^~j~ 1) = T(a+k)/T(a) denotes the Pochhammer symbol (Pochhammer [a, k] ) 

j=i 

or shifted factorial. is a hypergeometric term and fulfils the recurrence equation (k = 
0,1,...) 

_ (k + a 1 )...{k + a p ) 
k+1 - (k + b 1 )...(k + b q )(k + l) k 

with the initial value 

A := 1 . 

In Mathematica the generalized hypergeometric function is denoted by 
HypergeometricPFQ [plist , qlist , x] , where 

plist = [ai,a 2 , ■ ■ ■ ,a P ] and qlist = [&i, b 2 , b q ] . 

To find the representing Laurent-Puiseux expansion of a given expression / with respect to the 
variable x and point of development xq one uses the function call PowerSeries [f , {x,xO}] 
(short PS [f , {x,xO}] ). The package also supports the syntax PowerSeries [f ,x,xO] . If 
xq = one can also use either PowerSeries [f ,x] , or PS [f ,x]. Furthermore the procedure 
FunctionToHypergeometric [f ,x] brings the resulting series in hypergeometric notation. One 
has e.g. 

In [13]:= PS[ArcSin[x]~2,x] 

k 2 + 2 k 2 
4 x k! 

0ut[13]= sum[ , {k, 0, Infinity}] 

(1 + k) (1 + 2 k) ! 

In [14] := FunctionToHypergeometric [ArcSin [x] ,x] 

113 2 

Out [14]= x hypergeometricPFQ[{-, -}, {-}, x ] 

2 2 2 

In [15] := FunctionToHypergeometric [ArcSin [x] "2 ,x] 
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2 3 2 

Out [15]= x hypergeometricPFQ [{1, 1, 1>, {2, -}, x ] 

2 



In [16]:= PS[LaguerreL[n,x] ,x] 
k 

x Pochhammer [-n, k] 

0ut[16]= sum[ , {k, 0, Infinity}] 

2 

k! 

Here LaguerreL [n,x] denotes the n th Laguerre polynomial L n (x). 

We use the head sum to denote a resulting infinite sum since Mathematica Version 3 evaluates 
sums with head Sum instantly. Similarly Mathematica's builtin procedure HypergeometricPFQ 
has instant evaluations in some instances, e.g. 

In [17] := HypergeometricPFQ [{1,1}, {2}, x] 

Log[l - x] 
Out [17]= -( ) 

x 

By this reason we use the (non-evaluating) name hypergeometricPFQ instead. 

Note that the procedure PowerSeries can handle the more general case of Laurent-Puiseux 

series representations (p = 1,2,...) 



f( x ) = J2 a * 



k/p 



k=—oc 

for example 

In [18]:= PS[Sin[Sqrt[x]] ,x] 

k 1/2 + k 
(-1) x 

Out [18]= sum[ , {k, 0, Infinity}] 

(1 + 2 k) ! 



3 Recurrence Equations 

Similarly to the sum and product algorithms for holonomic differential equations there are 
simple sum and product algorithms for holonomic recurrence equations. These compute the 
recurrence equation for the sum and product, respectively, of two discrete functions and 
6^, given by two holonomic recurrence equations RE1 and RE2. These algorithms are invoked 
by the procedures SumRE [RE1 , RE2 , a [k] ] and ProductRE [RE1 ,RE2 ,a[k] ] , respectively. We 
get for example 
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In [19] := SumRE[a[k+l]==a[k]/(k+l) ,a[k+l] ==a[k] ,a[k]] 

2 

Out [19]= (1 + k) a[k] + (-1 - 3 k - k ) a[l + k] + 

> k (2 + k) a [2 + k] == 

In [20] := ProductRE[a[k+l]==a[k]/(k+l) ,a[k+l]==a[k] ,a[k]] 
Out [20]= a[k] + (-1 - k) a[l + k] == 

Here the first computation finds the recurrence equation valid for the sum of l/k\ (satifying 
a/t+i = a,k/(k + 1)) and 1 (satifying a/t+i = a^), and the second one for the corresponding 
product. Prom this point of view the result of the second computation is obvious. 
The automatic computation of holonomic recurrence equations for expressions is supported 
through the procedure HolonomicRE[expr,a[k]] . Examples are 

In [21] := HolonomicRE[(n!+k!"2)/k,a[n]] 

2 2 
Out [21]= (1 + n) a[n] + (-1 - 3 n - n ) a[l + n] + 

> n a [2 + n] == 

In [22] := HolonomicRE[(n!+k! ~2)/k,a[k]] 
3 

Out [22]= k (1 + k) (3 + k) a[k] - 

2 2 

> (1 + k) (1 + 3 k + k ) (3 + 3 k + k ) a[l + k] + 

2 

> k (2 + k) a[2 + k] ==0 

The work with special functions for which holonomic recurrence equations exist, is possible. 
In [23] := HolonomicRE[LegendreP[k,x] ,P[k]] 

Out [23]= (1 + k) P[k] - (3 + 2 k) x P[l + k] + 

> (2 + k) P[2 + k] == 

In [24] := HolonomicRE [Binomial [2*k,k] *HermiteH [k,x] ,a[k]] 
Out [24]= 8 (1 + 2 k) (3 + 2 k) a[k] 

> - 4 (3 + 2 k) x a[l + k] + (2 + k) a[2 + k] ==0 

One has also access to Zeilherger's algorithm. For this purpose the package uses an implemen- 
tation of Paule/Schorn from RISC [18]. Zeilberger's algorithm finds a holonomic recurrence 
equation for a sum 

S n := F ^ k ) 

k~k\ 
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if F(n, k) is a hypergeometric term w.r.t. both n and k, i.e., the term ratios 
F{n + l,k) ^ F{n,k + 1) . 

are both rational functions in n and k. One gets for example for the binomial powers 
In[25] := HolonomicRE [sum [Binomial [n,k] ~5,{k,0,n}-] ,a[n]] 

4 2 
Out [25]= 32 (1 + n) (292 + 253 n + 55 n ) a[n] + 

2 3 

> (-514048 - 1827064 n - 2682770 n - 2082073 n - 

4 5 6 

> 900543 n - 205799 n - 19415 n ) a[l + n] + 

2 3 

> (-79320 - 245586 n - 310827 n - 205949 n - 

4 5 6 

> 75498 n - 14553 n - 1155 n ) a[2 + n] + 

4 2 

> (3 + n) (94 + 143 n + 55 n ) a[3 + n] ==0 

Some time ago, such a result was worth a publication ([19], [4]). 
Similarly, for the Apery numbers 

one deduces the recurrence equation 

In [26] := HolonomicRE [sum [Binomial [n,k] ~2* 
Binomial [n+k , k] ~2 , {k , , n}] , A [n] ] 

3 

Out [26]= (1 + n) A[n] - 

2 

> (3 + 2 n) (39 + 51 n + 17 n ) A[l + n] + 

3 

> (2 + n) A [2 + n] == 

which played an essential role in Apery's proof (see e.g. van der Poorten's [20] entertaining 
presentation of Apery's proof and its history 1 ) of the irrationality of 

OO -I 

c(3) = £^- 

k=l 



1 Van der Poorten writes: "To convince ourselves of the validity of Apery's proof we need only complete the 
following exercise: show that the recurrence equation is valid for (1). Neither Cohen nor I had been able to 
prove this in the intervening two months. . . " 
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Zeilberger's algorithm can be used to prove identities. The Legendre polynomials P n (x), e.g. 
can be represented as any of 



P n (x) 



£ 

2Fl\ 



• n 
. k , 



-n — 1 \ / 1 — x 

k , 



-n, n + 1 



2n h w 



-n, — n 
1 



1 +a; 



x L»/2J 

on ' 



1 ~ x , 
2n - 2k\ 



,n-2k 



k=0 



t:) 2 p i 



k I \ n J 

n/2,-n/2 + 1/2 
-n + 1/2 



-n/2,-n/2 + 1/2 



a; 2 



By the following computations three of these representations are proved to be consistent 

In [27] := HolonomicRE [sum [Binomial [n,k] * 

Binomial [-n-1 , k] * ( ( 1-x) /2) ~k , {k , , n}] , P [n] ] 

Out [27]= (-1 - n) P[n] + (3 + 2 n) x P[l + n] + 

> (-2 - n) P[2 + n] ==0 

In[28] := HolonomicRE [l/2~n*sum [Binomial [n,k] ~2* 
(x-1) " (n-k) * (x+1) ~k , {k, ,n}] , P [n] ] 

Out [28]= (1 + n) P[n] - (3 + 2 n) x P[l + n] + 

> (2 + n) P[2 + n] == 

In[29] := HolonomicRE [l/2~n*sum[(-l) ~k*Binomial [n,k] * 

Binomial [2n-2k , n] *x~ (n-2k) , {k , , Infinity}] , P [n] ] 

Out [29]= (-1 - n) P[n] + (3 + 2 n) x P[l + n] + 

> (-2 - n) P[2 + n] ==0 

To finish the proof that these sums represent the same family of functions, it is enough to verify 
that two initial values agree since by the above computations the sums satisfy the same second 
order recurrence equation. Note that the last sum is supported in the interval [0, L^/2j] which 
can be seen from the two upper indices — n -^- and — ^ of the hypergeometric representation 
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In [30] := SumToHypergeometric[sum[l/2~n*(-l)~k*Binomial[n,k]* 
Binomial [2n-2k , n] *x~ (n-2k) , {k , , Infinity}] ] 

n 

Out [30]= (x Gamma [1 + 2 n] 

1-n-n 1 - 2 n -2 
> hypergeometricPFQ [{ , — }, { }, x ])\ 

2 2 2 



n 2 
> / (2 Gamma [1 + n] ) 

The procedure SumToHypergeometric [sum [expr , {k,0, Infinity}] ] converts an infinite sum 
into hypergeometric notation; the procedure FunctionToHypergeometric which was men- 
tioned before uses SumToHypergeometric internally. 

4 Generating functions 

Generating functions can be calculated by converting the holonomic recurrence equation of 
the coefficient sequence to the corresponding differential equation for its generating function 
if applicable. 

The procedure Convert [sum [expr , {k,k0, Infinity}] ,x] converts the Laurent-Puiseux se- 
ries Sum [expr, {k,k0, Infinity}] w.r.t. x to an expression if the corresponding differential 
equation can be solved by DSolve. For example, we get 

In[31]:= Convert [sum [Binomial [n,k] x~k,{k,0,lnf inity}] ,x] 

n 

Out [31]= (1 + x) 

In [32] := Convert [sum[x~ (2k+l) / (2k+l) ! ,{k,0, Infinity}] ,x] 

x 

-1 E 

Out [32]= + -- 

x 2 

2 E 

The procedure GeneratingFunction[a,k,z] tries to find the generating function 

oo 

Yj a k zk ) 
A;=0 

and the procedure Exponent ialGeneratingFunct ion [a, k,z] searches for the exponential 
generating function 

oo 

Both functions use Convert internally. We get e.g. for the Legendre polynomials 



9 



In [33]:= specf imprint 



In [34] := GeneratingFunction[LegendreP [n,x] ,n,z] 
SpecialFunctions , (C) Wolfram Koepf , version 1.00, 
> November 20, 1996 
specf un-info : RE: 

(1 + k)*a[k] - (3 + 2*k)*x*a[l + k] + 

(2 + k)*a[2 + k] ==0 
specf un-info : DE: 

2 

(-x + z) F[z] + (1 - 2 x z + z ) F' [z] == 
specf un-info : Trying to solve DE . . . 
specf un-info : DSolve computes 

C[l]/(1 - 2*x*z + z~2)~(l/2) 
specf un-info : expression rearranged: 

C[l]/(1 - 2*x*z + z~2)~(l/2) 
specf un-info : Calculation of initial values... 
specf un-info : C[l] == 1 



Out [34] = 



Sqrt[l - 2 x z + z ] 



The command specf unprint turns on a verbatim mode with which you will be informed 
about intermediate computation steps. 
Similarly, we get 

In[35] := ExponentialGeneratingFunction[HermiteH[n,x] ,n,z] 

SpecialFunctions, (C) Wolfram Koepf, version 1.00, 

> November 20, 1996 

specf un-info : REProduct entered 

specf un-info : RE: 

2*a[k] - 2*x*a[l + k] + (2 + k)*a[2 + k] ==0 
specf un-info : DE: 

2 (-x + z) f [z] + f ' [z] == 
specf un-info : Trying to solve DE . . . 
specf un-info : DSolve computes 

E~(2*x*z - z~2)*C[l] 
specf un-info : expression rearranged: 

E~((2*x - z)*z)*C[l] 
specf un-info : Calculation of initial values... 
specf un-info : C[l] == 1 



(2 x - z) z 

Out [35]= E 

In [36]:= nospecf unprint 



Note that nospecf unprint turns off the verbatim mode. 

You saw that Convert uses the builtin Mathematica procedure DSolve to solve the differential 
equation corresponding to the given series coefficients. This can be rather time consuming or 
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might be without success. If you know the generating function in advance, then you don't 

have to solve a differential equation. 

The following two statements prove, for example, that 



Y PJx) z n = , 



— k 

a k z 



modulo two initial values without using DSolve. Here the procedure SimpleRE[f ,x,a[k]] 
finds the recurrence equation which is valid for the series coefficients of f (see e.g. [12]) 

In [37] := SimpleRE [1/Sqrt [l-2*x*z+z~2] ,z,P[n]] 

Out [37]= -((1 + n) P[n]) + (3 + 2 n) x P[l + n] - 

> (2 + n) P[2 + n] == 

In [38] := HolonomicRE[LegendreP[n,x] ,P[n]] 

Out [38]= (1 + n) P[n] - (3 + 2 n) x P[l + n] + 

> (2 + n) P[2 + n] == 



5 Application: The Z- Transform 

The Z-Transform 

oo 

Z{a k } = f(z) = Y, 

of a sequence {a^} is the discrete analogue of the Laplace Transform. This series converges in 
the region outside the circle \z\ = \zq \ = limsup \/\ak\ ■ 

k—too 

The procedure ZTransf orm[a,k,z] computes the Z- Transform of the sequence a k , whereas 
the procedure InverseZTransf orm[f ,z] gives the inverse Z- Transform of the functions f(z), 
i.e., it computes a k . Internally these procedures use Convert and PowerSeries, respectively. 
As an application we give the examples 

In [39] := fun=ZTransf orm[l/k! ,k,z] 

1/z 

Out [39]= E 

In[40] := InverseZTransf orm[fun,z] 

1 k 
(-) 
z 

Out [40]= sum[ , {k, 0, Infinity}] 

k! 

In [41] := fun=ZTransform[l/(2*k+l) ! ,k,z] 
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Sqrt [1/z] 

-1 E 

Out [41]= + 

Sqrt [1/z] 1 1 

2 E Sqrt[-] 2 Sqrt [-] 

z z 

In[42] := InverseZTransf orm[fun,z] 

1 k 
(-) 
z 

Out[42]= sum[ , {k, 0, Infinity}] 

(1 + 2 k) ! 

In [43] := InverseZTransf orm[z/(l+z+z~2) ,z] 

1 k 1 1 + k 

(-( )) (-) 

-2 z 

1 + z 

Out [43]= sum[ , {k, 0, Infinity}] 

-2 

1 + z 



6 Application: Feynman Diagrams 

In [5], Fleischer and Tarasov gave the representation 2 

v(ty r , ) = ( _. )a+ p +1 r( a +/3+ 7 -d/2)r(d/2- 7 )r(«+ 7 -d/2)r(/3+ 7 -d/2) 

1 ' /,7; 1 ; r(a)r(/3)r(d/2)r(a + /3 + 2 7 -d)(m 2 ) a +P+"f- d 

fa + P + j-d, a + 7 -d/2 

■ 2-^1 I 

\ a + /3 + 2j -d 

for the calculation of a certain Feynman diagram. The function V(a, /3, 7 ) can be easily 
computed for a, (3, 7 G {0, 1}, hence one would like to attribute the calculation for nonnegative 
integer arguments a, (3, 7 to these cases. The simplest way is to compute the function by a 
pure recurrence equation that increases only one of its arguments, e.g. /3, and leaves the other 
arguments constant. By Zeilberger's algorithm, the package gives the pure recurrence equation 

In [44] := HolonomicRE[(-l)~(alpha+beta+gamma)* 

Gamma [alpha+beta+gamma-d/2] *Gamma [d/2-gamma] * 
Gamma [alpha+gamma-d/2] *Gamma [beta+gamma-d/2] / 
(Gamma [alpha] *Gamma [beta] *Gamma[d/2] * 
Gamma [alpha+beta+2*gamma-d] * (m~2) ~ 
(alpha+beta+gamma-d) ) * 
Hypergeometric2Fl [alpha+beta+gamma-d, 
alpha+gamma-d/ 2 , alpha+beta+2*gamma-d , z] , V [beta] ] 



2 I would like to thank Jochem Fleischer who informed me about a typographical error in formula (31) of [5]. 
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Out [44]= (2 beta - d + 2 gamma) 



> (2 alpha + 2 beta - d + 2 gamma) 

> (2+2 alpha + 2 beta - d + 2 gamma) V[beta] + 

2 

> 2 beta (2+2 alpha + 2 beta - d + 2 gamma) m 

> (2 alpha + 2 beta - 2 d + 4 gamma + 2 z + 

> 2 beta z - d z) V[l + beta] + 

> 8 beta (1 + beta) (1 + alpha + beta - d + gamma) 

4 

> m z V[2 + beta] == 

with respect to (3, and similar recurrence equations can be computed w.r.t. a and 7. 

7 Application: Polynomial Approximations of Solutions of 
Differential Equations 

Series solutions of holonomic differential equations satisfy holonomic recurrence equations from 
which their coefficients efficiently can be calculated by iteration. This is the asymptotically 
fastest known algorithm for the given purpose. It can be invoked by the Mathematica proce- 
dure SeriesSolution[DE,y [x] , initial value s ,n] . Here DE is a given differential equation 
in terms of y [x] , initialvalues is a list of initial values, and n is the order of the proposed 
approximation. For the Airy function we have for example 

In [45] := Timing [SeriesSolution[y' ' [x]-x*y[x]==0,y[x] , 
{AiryAi [0] , AiryAi ' [0] } , 10] ] 

x 

Out [45]= {0.3 Second, -( ) - 

1/3 1 
3 Gamma [-] 
3 

4 7 

x x 

> 

1/3 1 1/3 1 

12 3 Gamma [-] 504 3 Gamma [-] 
3 3 

10 

x 1 

> + + 

1/3 1 2/3 2 

45360 3 Gamma [-] 3 Gamma [-] 
3 3 
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3 6 

x x 

> + + 

2/3 2 2/3 2 

6 3 Gamma [-] 180 3 Gamma [-] 
3 3 

9 

x 

> } 

2/3 2 
12960 3 Gamma [-] 
3 

whereas the buitin Series command gives 

In [46] : = Timing [Series [AiryAi [x] , {x , , 10}] ] 

1 x 

Out [46]= {3.5 Second, + 

2/3 2 1/3 1 

3 Gamma [-] 3 Gamma [-] 
3 3 

3 4 
x x 

> + 

2/3 2 1/3 1 

6 3 Gamma [-] 12 3 Gamma [-] 
3 3 

6 7 
x x 

> + 

2/3 2 1/3 1 

180 3 Gamma [-] 504 3 Gamma [-] 
3 3 

9 10 

x x 

> _ + 

2/3 2 1/3 1 

12960 3 Gamma [-] 45360 3 Gamma [-] 
3 3 

11 

> [x] } 

Hence even for modest order the method can be rather fast. If the order is large, then 
the speedup is even more impressive as can be seen from the computations below. Here 
Taylor [f ,x,x0,n] gives a Taylor approximation of order n for f(x) at the point of develop- 
ment xq, by generating a holonomic differential equation for f(x), and using SeriesSolution, 
if applicable. 
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In [47] : = Timing [Series [AiryAi [x] , {x , , 100}] ; ] 
Out [47]= {33.1833 Second, Null} 
In [48] : = Timing [Taylor [AiryAi [x] , -fx , , 100}] ; ] 
Out [48]= {2.03333 Second, Null} 

In [49] : = Timing [Series [Sin [x] *Exp [x] , {x , , 100}] ; ] 
Out [49]= {5.88333 Second, Null} 

In [50] : = Timing [Taylor [Sin [x] *Exp [x] , {x , , 100}] ; ] 
Out [50]= {1.11667 Second, Null} 

In [51] := Timing [Series [Sin [x]*Exp[x] ,{x, 0,1000}] ;] 
Out [51]= {8517.53 Second, Null} 

In [52] : = Timing [Taylor [Sin [x] *Exp [x] , {x , , 1000}] ; ] 

Out [52]= {7.63333 Second, Null} 

In [53]:= ps=PS [Sin[x] *Exp [x] ,x] 

k/2 k k Pi 

2 x Sin[ ] 

4 

Out [53]= sum[ , {k, 0, Infinity}] 

k! 

In[54] := Timing[ps/ . {sum->Sum, Inf inity->1000} ;] 
Out [54]= {7.21667 Second, Null} 

Note that even the closed form solution given in line 43 does not provide a faster way to 
compute the series: The evaluation of each single coefficient takes about the same time than 
the iterative computation of the coefficients used by Taylor. 

8 Advanced Applications 

In this section, we give some advanced applications. Without explicit mentioning, all compu- 
tations give proofs modulo initial values. 
The first one is DougaWs identity 



7-^6 



a, 1 + |, b, c, d, 1 + 2a — b — c — d + n, —n 
|, 1 + a — b, 1 + a — c, 1 + a — d, b+c+d—a — n, 1 + a+n 



1 



(I + a) n (l + a - b - c) n {\ + a - b - d) n {\ + a - c - d) n 
(1 + a - + a - c) n (1 + a - d) n (1 + a - b - c - d) n 

15 



which is proved by 



In [55] : 



HolonomicRE [sum[HyperTerm[ 

{a, l+a/2,b,c,d, l+2*a-b-c-d+n, -n} , 

{a/2 , 1+a-b , 1+a-c , 1+a-d , b+c+d-a-n , 1+a+n} , 1 , k] , 

{k , -Inf inity , Inf inity}] , S [n] ] 



Out [55] 



(1 + a + n) (1 + a- b- c + n) 



> 



(1 + a 



b - d + n) (1 + a - c - d + n) S[n] - 



> 



(1 + a 



b + n) (1 + a - c + n) (1 + a - d + n) 



> 



(1 + a 



b - c - d + n) S[l + n] == 



Note that from this computation the parameters of the right hand hypergeometric term in (2) 
can be directly read off. Hence Zeilberger's algorithm has discovered the right hand side of 
(2) from its left hand input. Here the Mathematica procedure HyperTerm[plist ,qlist ,x,k] 
denotes the kth summand of the series HypergeometricPFQ [plist ,qlist ,x] . 
Similarly Clausen 's identity- 



is proved by the computation 

In [56] := HolonomicRE [HypergeometricPFQ [ 
{a,b, 1/2-a-b-n, -n} , 
{l/2+a+b,l-a-n,l-b-n},l] ,S[n]] 

Out [56]= -((2 a + n) (a + b + n) (2 b + n) S [n] ) + 

> (a + n) (b + n) (2 a + 2 b + n) S[l + n] == 

The only case when the square of a 2-F1 function is a 3F2 function is given by Clausen's formula 



Rather than proving this in the usual way by showing that both sides satisfy the same holo- 
nomic differential equation, we can use Zeilberger's algorihm to determine the right hand 
side of (3), given the left hand product, by using the Cauchy product. This is done by the 
computation 

In[57]:= HolonomicRE [sum [ 




a, 6, 1/2 — a — b — n, —n 
1/2 + a + 6, 1 — a — n, 1 — b — n 

= (2a) n (a + b) n (2b) n 
{2a + 2b) n {a) n {b) n 





HyperTerm [{ a , b} , {a+b+ 1/2} , 1 , j ] * 
HyperTerm[{a,b},{a+b+l/2}, l,k-j] , 
{j,0,k>] ,a[k]] 
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Out [57]= -2 (2 a + k) (a + b + k) (2 b + k) a[k] + 



> (l + k)(2a + 2b + k)(l + 2a + 2b + 2k) 

> a[l + k] == 

Observe that one can directly read off the parameters of the hypergeometric function on the 
right hand side of (3). 

The following identity of Askey and Gasper 

(a + 2)„ ( -n, n + a + 2, 2±i ' 



n! 



3-^2 



[ ^ ] (ti. ? (f + 1 )»- J (^)„- 2? ("+ 1 )^ 



2j-n,n-2j + a + l,^ 
a + l,a±2 



a+l 

.< Z 7 — //,. // — Z 7 -t- f r — t— I . 

(see e.g. [7]) was an essential tool in the proof of the Biebeibach conjecture by de Branges [3]. 
The Askey-Gasper identity is proved by the computation 

In [58] := HolonomicRE [sum [Pochhammer [1/2 , j] * 

Pochhammer [alpha/2+1 ,n-j] *Pochhammer [(alpha+3) II ,n-2* j] * 
Pochhammer [alpha+1 ,n-2* j] /Pochhammer [(alpha+3) /2 ,n-j] / 
Pochhammer[(alpha+l)/2,n-2*j]/(n-2*j) ! / j ! * 
HyperTerm[{2*j -n ,n-2* j+alpha+1 , (alpha+1) /2}- , 
{alpha+1, (alpha+2)/2},l,k] ,{j , -Infinity, Infinity}] ,a[k]] 

Out [58]= (1 + alpha + 2 k) (k - n) (2 + alpha + k + n) 

> a[k] - (1 + k) (1 + alpha + k) 

> (3 + alpha + 2 k) a[l + k] ==0 

Note that here the input is much more complicated than the output is! The computational 
trick is that we interchanged the order of summation of the double sum. 
By an application of Clausen's formula from the Askey-Gasper identity it follows that for 
a > — 2 the left hand function in (4) is nonnegative. This was the result which de Branges 
needed for a = 0, 2, . . . for his proof of the Bieberbach conjecture. 

In recent work the parameter derivatives for the Jacobi, Gegenbauer and Laguerre polynomials 
have been determined ([6], [14]). For the Laguerre polynomials lI?\x) one gets for example 
the simple formula 

71-1 



k=0 

This result can be obtained, e.g., by taking \i — > in the connection relation 



The latter statement can be discovered by 
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In [59] := HolonomicRE [sum[Pochhammer [mu,n-k] / (n-k) !* 
(-1) ! *Binomial [k+alpha,k-j] *x~j , 
{k , -Inf inity , Inf inity}] , a [ j ] ] 

Out [59]= (-j + n) x a[j] + (1 + j) 

> (1 + alpha + j + mu) a[l + j] ==0 

How can we obtain the exponential generating function 

£ - P n {x) z n = e xz Jo (zVTl^) (5) 

n=0 

of the Legendre polynomials, J n (x) denoting the Bessel functions? Note that Jo(x) is a o-^i 
hypergeometric function 

In [60] := FunctionToHypergeometric[BesselJ[0,x] ,x] 

2 

-x 

Out [60]= hypergeometricPFQ [{}, {1}, ] 

4 

To obtain (5), we use the representation 

1 \ 



/ -n/2, (1 -n)/2 
P n (x)=x n 2 FA 1 " 



1 2 

x z 



for the Legendre polynomials, and change the order of summation: 

In[61]:= HolonomicRE [sum [x~n* 

HyperTerm [{-n/2 , ( 1-n) /2} , { 1} , 1- l/x~2 , k] /n ! *z~n , 
-[n , -Inf inity , Inf inity}] , a [k] ] 

2 2 
Out [61]= (1 - x) x (1 + x) z a[k] + 

2 2 

> 4 (1 + k) x a[l + k] == 



Now we can read off the hypergeometric parameters of the Bessel function, and using the 
initial value 



n=0 



we obtain (5). 
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9 Efficient Computation of Orthogonal Polynomials 



In recent studies we discovered that none of the popular computer algebra systems Axiom, 
Macsyma, Maple, Mathematica, Mupad and REDUCE had implemented the most efficient 
algorithms to compute specific orthogonal polynomials, neither for rational exact nor for 
numerical computations, although Mathematica came rather near [16]. 

Hence, in the package SpecialFunctions, we have included efficient implementations for this 
purpose. If you compute the 1,000th Chebyshev polynomial by Timing [ChebyshevT [1000 , x] ; ] 
before and after loading the package, you will see how the computation is accelerated using 
the package. The larger n is, the larger the advantage. (If you have a fast computer and 
enough memory, you should try n = 5, 000 or n = 10, 000.) 

For n = 0, 1, 2, . . . the package supports the computation of the classical orthogonal polynomi- 
als JacobiP [n,a,/3,x] , GegenbauerC [n,a,x] , LaguerreL[n,a,x] and HermiteH[n,x] , as 
well as the computation of the discrete orthogonal polynomials (see [17]) Hahn[n,N,a,/3,x] , 
DiscreteChebyshev[n,N,x] , Meixner [n,7,/L/,x] , Krawtchouk[n,N,p,x] , Discrete- 
Laguerre [n , p , a , x] and Charlier [n,/u ,x] . Note that these as well as more special func- 
tions can also be used in connection with HolonomicDE, HolonomicRE, etc. 
With ?SpecialFunctions ' * you get a list of all procedures exported by the package, and by 
?PowerSeries, e.g., you get a description of the Power Series procedure, together with an 
example. Moreover the contents of the integrated package fastZeil by Paule/Schorn can be 
listed by TfastZeil'*. 

Some examples for the computation of specific orthogonal polynomials: 

In[62]:= Charlier [5,mu,x] 

5 x 10 (1 - x) x 

Out [62]= 1 

mu 2 
mu 

10 (1 - x) (2 - x) x 

> _ 

3 

mu 

5 (1 - x) (2 - x) (3 - x) x 

> 

4 

mu 

(1 - x) (2 - x) (3 - x) (4 - x) x 

> 

5 

mu 

In [63]:= ChebyshevT [1000 , 1/4] 
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Out [63] = 463388825262072952755531624295890954893293543\ 

> 931057925304868803092234816485137547165153207603\ 

> 480107997324539642132129419798682743029733063451\ 

> 024106656535321705939115104797434131268364882269\ 

> 510117040295963467578579678911640493759116734331\ 

> 488874299433540209147382435657218918439668245092\ 

> 8217512771529377 / 

> 2143017214372534641896850098120003621122809623411\ 

> 067214887500776740702102249872244986396757631391\ 

> 716255189345835106293650374290571384628087196915\ 

> 514939714960786913554964846197084214921012474228\ 

> 375590836430609294996716388253479753511833108789\ 

> 215412582914239295537308433532085966330524877367\ 

> 4411336138752 

In [64] := N[ChebyshevT[10~10,N[l/4, 100]] ,100] 

Out [64] = -0 . 161590710058830973064545131784268650111183\ 

> 53081684739911777110975137352756666312533037757 

Note that for plotting purposes, Mathematical original implementation is advantageous. 
Using the Charlier polynomials, we finally give some examples of the symbolic use of 
discrete orthogonal families: 

In[65] := HolonomicRE [Charlier [n,mu,x] ,C[x]] 

Out [65]= (-1 - x) C[x] + (1 + mu - n + x) C[l + x] - 

> mu C[2 + x] == 

In [66] := HolonomicRE [Charlier [n,mu,x] ,C[n]] 

Out [66]= (1 + n) C[n] + (-1 - mu - n + x) C[l + n] + 

> mu C[2 + n] == 

In [67] := HolonomicRE [Binomial [mu,n] *Charlier [n,mu,x] ,a[n]] 
Out [67]= (-mu + n) (1 - mu + n) a[n] + 
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> (1 - mu + n) (1 + mu + n - x) a[l + n] + 

> mu (2 + n) a [2 + n] == 

In [68] := HolonomicRE [Charlier [n,mu,2x] ,a[x]] 
Out [68]= 2 (-3 - mu + n - 2 x) (1 + x) (1 + 2 x) 

2 3 

> a[x] +(6+2mu+mu +mu -lln-7mun- 

2 2 2 3 

> 3mun + 6n+3mun-ii+22x + 

2 

> 6mux + 2mu x-24nx-8munx+ 

2 2 2 2 3 

> 6nx+24x+4mux-12nx+8x) 

2 

> a[l + x] + mu (-1 - mu + n - 2 x) a [2 + x] == 
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